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' Results from hydrodynamical SPH simulations of galaxy clusters are used to 

^ ' investigate the dependence of the final cluster X-ray properties on the numerical 

^ _ resolution and the assumed models for the physical gas processes. Two differ- 

Xf^ • ent spatially fiat cosmological models have been considered: a low-density cold 

^ ■ dark matter universe with a vacuum energy density = 0.7 (ACDM) and a 

O ' cold+hot dark matter model (CHDM). For each of these models two different 

clusters have been extracted from a cosmological A^— body simulation. A series 
of hydrodynamical simulations has then been performed for each of them using a 
■ TREESPH code. These simulations first include radiative cooling and then also 

I ' conversion of cold gas particles into stars; because of supernova explosions these 

^ • particles can release energy in the form of thermal energy to the surrounding 

^ ■ intracluster gas. For a specific treatment for the thermal state of the gas, simu- 

lation runs have been performed with different numerical resolutions. This is in 
order to disentangle in the final results for the cluster profiles, the effects of the 
j_| ■ resolution from those due to the assumed model for the gas thermal evolution. 

^ The numerical resolution of the simulation is controlled by the number of gas par- 

ticles Ng and the chosen value for the gas gravitational softening parameter Eg. 
The latter is proportional to the minimum SPH smoothing length and therefore 
sets a maximum spatial resolution for the simulations. For the cooling runs, final 
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X-ray luminosities have been found to be diverging according to Lx oc 1/e 
The gas density profiles are also diverging at the cluster center. This is in agree- 
ment with previous findings. When cold gas particles are allowed to convert into 
stars, the divergences are removed. The final gas profiles show a well defined 
core radius and the temperature profiles are nearly flat. For the most massive 
test cluster in the ACDM model, these simulations show a prominent cooling 
flow in the cluster core. This cluster was analyzed in detail, running simulations 
with different star formation methods and increasing numerical resolution. A 
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comparison between different runs shows that the results of simulations, based 
on star formation methods in which gas conversion into stars is controlled by 
an efficiency parameter c^, are sensitive to the numerical resolution of the simu- 
lation. In this respect star formation methods based instead on a local density 
threshold, as in Navarro and White (1993), are shown to give more stable results. 
Final X-ray luminosities are found to be numerically stable, with uncertainties 
of a factor ~ 2. These simulations are also in good agreement with observational 
data when the final results are compared with the observed star formation rate 
and the luminosity-temperature relation from cooling flow clusters. Therefore I 
find that hydrodynamical simulations of cooling clusters can be used to give reli- 
ably predictions of the cluster X-ray properties. For a given numerical resolution, 
the conversion of cool gas particles into stars as in Navarro and White should be 
preferred. 

Subject headings: cosmology: clusters-galaxies: clusters-methods:numerical 

1. Introduction 

Galaxy clusters are the largest virialized structures known in the universe. According 
to the hierarchical scenario their evolution rate is a strong function of the background cos- 
mology (Peebles, Daly & Juszkiewicz 1989; Lilje 1992; Viana & Liddlc 1999; Oukbir & 
Blanchard 1992; Eke et al. 1998; Sadat, Blanchard & Oukbir 1998; BahcaU, Fan & Cen 
1997), thus making galaxy clusters natural tools for constraining the cosmological models. 
Galaxy clusters are also powerful X-ray emitters. X-ray observations have shown that most 
of the baryons in galaxy clusters are in the form of hot ( T ~ 10^ ) ionized X-ray emit- 
ting gas (Forman & Jones 1982; Sarazin 1986). The bulk of the emission is via thermal 
bremsstrahlung; its dependence on the square of the gas density allows one to select cluster 
samples without the contamination effects which may arise in the optical band. For this 
reason, galaxy clusters have been the subject of extensive observational programs in the 
X-ray band (Henry et al. 1992; Henry 1997; Ebeling et al. 1997; Rosati et al. 1998). 
Observations of the cluster X-ray emission and temperature can be used to reconstruct the 
radial gas density and temperature profile, assuming spherical symmetry. The gas density 
profile has a radial fall-off with slope ~ —2 and a constant density core in the inner regions, 
with typical size Vg ~ 50 — 200h~^Kpc (Sarazin 1986); the gas temperature profile is instead 
nearly isothermal within the virial radius. The density of the gas can be used to recover 
the dark matter profile and the total cluster mass. Assuming virial equilibrium these clus- 
ter properties are connected to the primeval power spectrum and the assumed cosmological 
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model, thus providing important clues for testing cosmologies (Makino, Sasaki & Suto 1998; 
Navarro, Prenk & White 1997). 

Other cosmological information can be obtained from the statistical properties of the 
ensemble of X-ray clusters. X-ray observations of cluster number counts, the X-ray tem- 
perature function (Henry & Arnoud 1991; Edge et al. 1990; Henry 1997) and the X-ray 
luminosity function (Rosati et al. 1998; Ebeling ct al. 1998) arc powerful probes for 
constraining the values of the cosmological parameters Qq and (Xg (Henry & Arnoud 1991; 
White, Efstathiou & Frenk 1993; Bahcall & Fan 1998; Eke, Cole & Frenk 1996; Kitayama, 
Sasaki & Suto 1998). An analytical framework for connecting the X-ray temperature and 
the luminosity function to theoretical models can be obtained, within the Press & Schecter 
(1974) approach, by assuming hydrostatic equilibrium for the gas distribution and neglecting 
radiative cooling. In accordance with these assumptions, analytical methods can then be 
used to derive predictions for the evolution of the X-ray luminosity function and its correla- 
tion with the cluster X-ray temperature (Kaiser 1986; Kitayama & Suto 1996; Mathiesen 
& Evrard 1998; Viana & Liddle 1999). 

Given the wealth of information which can be obtained from observations of X-ray 
clusters, a lot of efforts have been devoted to obtaining directly the gas and temperature 
distributions, using numerical simulations for investigating the evolution of galaxy clusters. 
Collisionless N-body simulations have been used to study substructure formation (West, 
Oemler & Dekel 1998; Crone, Evrard & Richstone 1996; Buote & Xu 1997; Jing et al. 
1995), density profiles (Crone, Evrard & Richstone 1994; Huss, Jain & Steinmetz 1999; 
Moore et al. 1998; Jing & Suto 2000), and statistical properties in a cosmological volume 
(Eke, Cole & Prenk 1996; Lacey & Cole 1994). In these simulations the dark matter 
distribution is a good tracer of that of the gas provided that the cluster dynamical state 
is not far from equilibrium, with a small fraction of substructures within the virial radius. 
Furthermore, the dark matter density profile shows no evidence for a core radius (Moore 
et al. 1998) whatever numerical resolution is achieved. Numerical simulations have then 
been extended to include the hydrodynamics of the gaseous component. The numerical 
methods used are either Eulerian (Cen & Ostriker 1994; Anninos & Norman 1996; Brian 
& Norman 1998; Kang et al. 1994; Bryan et al. 1994; Cen 1997; Cen et al. 1995), with 
a fixed or adaptative grid, or Lagrangian (Evrard 1988, 1990; Thomas & Couchman 1992; 
Navarro, Prenk & White 1995; Eke, Navarro & Prenk 1998; Katz & White 1993; Yoshikawa, 
Itoh & Suto 1998; Valdarnini, Ghizzardi & Bonometto 1999). The Lagrangian schemes 
are based on the smoothed particle hydrodynamic technique (SPH: Gingold & Monaghan 
1977; Monaghan 1992). In these simulations the required dynamical range can be quite 
demanding. A large simulation box is needed in order to obtain a meaningful statistical 
sample of galaxy clusters; at the same time the minimum spatial resolution must be at 
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least close to the cluster core radii, where the bulk of the X-ray emission originates. For 
Lagrangian methods a partial solution is the multimass technique (Katz & White 1993), 
where the particles increase their masses according to their distance from the cluster center. 
A single cluster can then be simulated with a comparatively high numerical resolution, with 
external shells of matter surrounding the cluster and representing the large-scale gravitational 
field. In these simulations the gas component is treated as a single adiabatic fluid, without 
taking into account the effects of radiative cooling, and the physical processes of merging, 
substructure formation, shocks and compressional heating of the gas can be modelled in this 
way. A comparison between different numerical simulations shows that they are successful 
in reproducing the gross features of the cluster properties (Prenk et al. 1999). 

The inclusion of radiative cooling for the gas is important on scales where the gas cooling 
time is shorter than the Hubble time. For galaxy clusters the spatial extent of this region 
is between 50 and 200 Kpc from the cluster center. A gas cooling radius TcooI can then be 
defined where the two time scales are equal. There is of observational evidence that numerical 
simulations must include radiative cooling, together with star formation and energy feedback, 
in order to model adequately the relevant physical processes of the gas during its evolution. 
In the inner regions where the radiative time scales are short, a cooling flow develops (Fabian 
1994), with Tcooi — 50h~^Kpc. More than 50% of clusters are estimated to have their cores 
in this phase (Peres et al. 1998). These instabilities can affect several properties, hke the 
Lx-Tx relation (Fabian et al. 1994; Allen & Fabian 1998a; Markevitch 1998). In addition 
to cooling flows there is also strong observational support for non-gravitational heating of the 
cluster gas. Simple scahng arguments predict Lx oc T|- (Kaiser 1986), while the observed 
relation satisfles Lx o: T|- (David et al. 1993). This discrepancy has been suggested 
by many authors (Evrard & Henry 1991; Wu, Fabian & Nulsen 1998; Ponman, Cannon & 
Navarro 1999) as being evidence for substantial heating of the intracluster gas due to energy 
injection by supernova explosions at high rcdshifts. Another central question for including 
star formation in simulations of the gaseous component is the growing observational evidence 
for radial gradients in the iron abundance (Ezawa et al. 1997; Fukazawa et al. 1998), with 
possible connections to coohng flows (Allen & Fabian 1998b). 

With increasing availability of computational power, numerical hydrodynamical simula- 
tions have attempted to model the effects of radiative coohng of the gas in the formation and 
evolution of cluster galaxies (Katz & White 1993; Suginohara & Ostriker 1998; Anninos & 
Norman 1996; Yoshikawa, Jing & Suto 2000; Pearce et al. 2000; Lewis et al. 2000). The 
numerical problems posed by the inclusion of gas cooling are challenging, mainly because the 
required increase in spatial resolution also requires that one keeps two-body heating mech- 
anisms under control. Moreover, it will be seen that the inclusion of radiative cooling for 
the gas cannot be separated from considering also star formation and energy feedback from 
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SN explosions, in order to obtain realistic cluster density profiles and luminosities. Previous 
simulations have produced some conflicting results (Yoshikawa, Jing & Suto 2000; Pearce 
et al. 2000), and so the question of the minimum resolution in this kind of simulation is 
still to be fully settled. The purpose of this paper is to test the numerical reliability of SPH 
hydrodynamical simulations of cluster formation. These simulations will include the effects 
of radiative losses, star formation and energy feedback from SN. Four different clusters from 
two cosmological models have been studied in several simulations with different numerical 
inputs and star formation modelling. Final profiles are compared in order to assess the ef- 
fects on the integrations of numerical resolution, and different star formation prescriptions. 
The paper is organized as follows. In §2 I describe the hydrodynamical simulations with 
radiative cooling and star formation that have been performed. The simulation results are 
then discussed in §3. In particular, §3.1 is dedicated to a comparison between different runs 
of the final radial density and temperature profiles, as well as of the X-ray luminosities. In 
§3.2 the simulation results are compared with previous findings. In §3.3 simulation runs with 
different star formation prescriptions are performed for a chosen test cluster which showed 
a well defined cooling instability. In §3.4 simulation results for runs with different prescrip- 
tions for star formation are compared against observed data from cooling flow clusters, in 
order to assess the consistency with real data for the assumed star formation models in the 
simulations. Finally the main conclusions are summarized in §4. 



2. Simulations 



In a previous paper (Valdarnini, Ghizzardi & Bonometto 1999, hereafter VGB) a large 
set of hydrodynamical simulations was used to study global X-ray cluster morphology and its 
evolution. The simulations were run using a TREESPH code with no gas cooling or heating. 
In order to assess the numerical reliability of the numerical integrations, four different clusters 
were selected as a representative sample of all of the simulation clusters. For this cluster 
sample, a large set of different integrations was performed by varying two numerical input 
parameters: the number of particles and the softening parameter. A comparison between the 
flnal gas density and temperature proflles, as well as with X-ray luminosities, then allowed 
the two-body heating to be kept under control and a fairly safe range of allowed values for the 
numerical parameters to be established. The simulation tests showed the relative importance 
of different numerical effects in these adiabatic simulations. It is therefore natural to use the 
same cluster sample to study the effects of including additional physics such as gas cooling 
and star formation. A comparison with the previous tests in VGB will show, with respect 
to the adiabatic case, the effects on the final cluster properties of considering energy sinks 
and non-gravitational heating. Here I give a short description of the cluster simulations 
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performed in VGB; the reader is referred to the original paper for more details. Numerical 
modelling of gas cooling and gas conversion into coUisionless stars is described later. 

In VGB three spatially flat cosmological models have been considered. A standard cold 
dark matter model (CDM) , a vacuum-energy dominated model with = 0.7 (ACDM) and 
a mixed dark matter model (CHDM). For the Hubble constant Hq — lOOhKm sec~^ Mpc"^, 
h = 0.5 is used for CDM and CHDM and h — 0.7 for ACDM. For all models, the primeval 
spectral index n = 1 and the baryon density parameter Q^h^ = 0.015. For CHDM = 0.20 
is the HDM density parameter of massive neutrinos; only one massive species is considered. 
All of the models were normalized in order to reproduce the present cluster abundance 
(Eke, Cole & Frenk 1996; Girardi et al. 1998). The cosmologies were chosen in order 
to have simple models with different properties. For each model an N-body cosmological 
simulation was first run in a L = 200h'~^Mpc comoving box using a P^M code. The particle 
number were Np = 10^ for the CDM and CHDM models with = 1, while Np = 84^ for 
ACDM with = 0.3. The same random numbers were used to set the initial conditions 
for all three cosmological models. The simulations started from an initial redshift Zm- At 
z = clusters of galaxies were located using a friends-of-friends (FoF) algorithm, so as to 
detect overdensities in excess of ~ 200Q^'^. For statistical analysis, VGB selected for each 
model the 40 most massive clusters. For each of these clusters a TREESPH hydrodynamical 
simulation was performed in physical coordinates. The integration was accomplished by 
first locating at z — the cluster center and identifying all of the simulation particles of 
the cosmological simulation within r2oo, where the cluster density is ~ 200Q~°-^ times the 
background density. These particles are located back at Zin, in the original simulation box, 
and a cube of size Lc enclosing all of them is then found, with a size ~ 15 — 25h~^ Mpc. A 
lattice of Nl = 22? grid points is set inside this cube; different lattices were used for each 
matter component. At each node position is associated a particle of corresponding mass and 
coordinates. The particles were then perturbed, using the same random realization as for 
the cosmological simulations. Additional frequencies are introduced so to sample the higher 
Nyquist frequency. The baryon particles are perturbed identically to the CDM particles 
and their initial temperature is set to % — 10^ °K. For the TREESPH simulations all the 
particles which lie inside a sphere of radius ivc/2 are kept. External gravitational fields are 
modelled by considering a larger cube of side 2Lc, inside the cube particle positions are set as 
for the smaller cube, but with no gas. The number of grid points is the same as for the inner 
cube, so that masses are 8 times larger than those of the inner cube; particles of the larger 
cube are considered only outside the smaller one. After the particle positions are perturbed, 
only those within a sphere of radius Lc from the cluster center are kept for the TREEPSH 
simulations. This multimass grid technique has already been used in cluster simulations by 
Katz & White (1993) and Navarro, Frenk and White (1995). 
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For each particle, the gravitational softening parameters are set according to the scaling 

1/3 

ffj oc m/ . The numerical integrations were performed with a tolerance parameter ^ = 0.7, 
without quadrupole corrections. The reliability of the numerical resolution was tested in 
VGB by taking the most massive and least massive clusters (labels 00 and 39, respectively) 
for the two models CHDM and ACDM. For each of them, different numerical tests were 
carried out. Accordingly to VGB, the cluster simulations can be performed with an ade- 
quate resolution using a number of gas particles Ng ^ 5, 000 and a gas softening parameter 
Eg ^ 50h~^Kpc. For the coUisionless component the corresponding values are scaled accord- 
ingly. In Table 1 the reference values for the four clusters are given. If one introduces 
radiative cooling then the final values for the core radius of the gas density are expected to 
be smaller than in the no-cooling case. This implies that the numerical simulations must 
have smaller values for the gas softening Eg, which in SPH sets a minimum resolvable scale, 
since SPH smoothing lengths are constrained to be ^ The shape of the density profile 

implies, with respect to the adiabatic case, higher values for the gas and dark matter den- 
sities at the core radius. This in turn implies that larger values for the number of particles 
are required in simulations of cooling clusters. This is essential in order to reduce the values 
of particle masses and hence the two-body heating time r^, which approximately scales as 
Tr oc l/{pdmd), where pd and are the dark matter density and particle mass. 

In order to check for these numerical effects, for each of the four test clusters a set of 
five TREESPH simulations was performed, using the same initial conditions and with the 
inclusion of radiative cooling, but for different values of Eg and Ng. With respect to the 
reference case, the other matter components have their particle numbers changed in propor- 
tion to the change in Ng. For the cluster ACDMOO, Table 2 reports the values of Eg and 
Ng for the five simulation runs. Table 4 is for CHDMOO. For these two clusters the generic 
simulation has cluster index clOO — j, with j = 00, 01, 05. The cluster clOO — 00 is the 
reference case without cooling. The simulations have been carried out with the same values 
for the other numerical parameters that were used in VGB, with the difference that here the 
minimum allowed time step for gas particles is Atm — 6.9 lO^yr. For the clusters ACDM39 
and CHDM39, the parameters of the numerical tests cl39 — j are given in Table 5. For these 
simulation runs, the number of particles and the initial redshifts are the same as for the 00 
clusters. Therefore Table 5 reports only the particle masses and softenings. The effects of 
radiative cooling are modelled in these simulations by adding to the SPH thermal energy 
equation an energy-sink term (Hernquist & Katz 1989, eq. [2.29]). The total cooling func- 
tion includes contributions from recombination and collisional excitation, bremsstrahlung 
and inverse Compton cooling. For cluster temperatures which satisfy T ^ 2KeV, a con- 
dition which is always satisfied for the cluster sample studied here, the dominant cooling 
mechanisms are free-free transitions, but line cooling becomes important for small clusters 
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and groups. For the same reason, heating from an ionizing UV background is not included in 

the thermal energy equation. The radiative cooling is computed for a gas having primordial 
abundances X = 0.75, Y = 0.25 with zero metallicity. In the simulation runs where the gas 
particles are allowed to convert part of their mass into stars (see below), the back effects 
of metallicities on the cooling function are neglected. This is a valid approximation as long 
as T ^ 2KeV and stellar metallicities are below Z ^ Zq {[Fe/H] ^ 0) (Brian & Norman 
1998; Carraro, Lia & Chiosi 1998). 

In addition to these simulations, which include the effects of radiative cooling, a mirror 
simulation was performed for each of them. The mirror runs had an additional prescription 
which allowed eligible gas to be turned into stars. These simulations are indexed as clOO — k 
and cl39 — k, where k — j + 5 and j is the index of the pure coohng runs in Tables 2, 4 & 5. 
For the cluster ACDMOO the simulation parameters for the cooling runs with star formation, 
clOO — k, are given in Table 3. The numerical parameters of the simulations with radiative 
cooling and star formation are the same as the corresponding cooling simulations. For this 
reason Table 3 reports parameters only for the cluster ACDMOO. The additional simulation 
with k = 11 in Table 3 has the same parameters as the A; = 10 run, but with 6 = 1 and 
quadrupole corrections, instead of ^ = 0.7. This is done in order to check the accuracy of the 
gravitational integration when a coUisionless population, with a different distribution from 
that of dark matter, is added to the simulation. Allowing the gas to cool radiatively will 
produce dense clumps of gas at low temperatures (~ 10^ °K). Thus cooling times will become 
shorter and even denser regions will develop. This is known as the overcooling instability 
(Suginohara & Ostriker 1998). In these regions the gas will be thermally unstable and will 
meet the physical conditions to form stars. In TREESPH simulations, star formation (SF) 
processes have been implemented using two different algorithms (Katz 1992; Navarro & 
White 1993). According to Katz (1992), a gas particle is in a star forming region if the flow 
is convergent and the local sound crossing time is larger than the dynamical time (i.e. the 
fluid is Jeans unstable). These two conditions read 



where Vi is the particle velocity, hi the SPH smoothing length and q is the local sound 
velocity. In a more refined version, Katz, Weinberg & Hernquist (1996, hereafter KWH) 
introduced two additional requirements: a star forming region must have a minimum physical 
hydrogen number density uh — O.lcm"^ and the local gas density must satisfy Pg/pg > 55.7 
( this follows from an isothermal profile giving a mean virialized over density of ~ 169). If 
a gas particle meets these criteria then it is selected as an eligible particle to form stars. 
In regions where the gas density is depressed because of gravitational softening, the Jeans 
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criterion is not applied, in order to avoid an underestimate of the local star formation rate 
(SFR) (Katz 1992, eq.[2]). The local SFR obeys the equation 

dpg/dt = -c^Pg/Tg = -dp^/dt , (2) 

where Pg is the gas density, p^ is the star density , Cy^ is a characteristic dimensionless efficiency 
parameter, Tg is the local collapse time (the maximum of the local cooling time and the 
dynamical time Td). Gas particles with T ^ 10^ °K have long cooling times and Tg — Td- The 
probability that a gas particle will form stars in a time step At is given by 

p = 1 - exp(-c^Ai/Tg). (3) 

A uniform random number is generated at every time step for each of the gas particles 
satisfying the star formation criterion, and equation (3) is used to compute the formation 
probability p. If < P then a mass fraction e^, of the mass of the gas is converted into a new 
coUisionless particle. This star particle has the position, velocity and gravitational softening 
of the original gas particle. Typical assumed values are — 1/3 and = 0.1 (KWH). 

The second algorithm for implementing SF in TREESPH simulations has been intro- 
duced by Navarro & White (1993 , hereafter NW). According to NW, any gas particle which 
is in a convergent flow and for which the density exceeds a threshold , i.e. 

( V-Vi< 

\ Pg > Pg,c = 7 ■ 10-2^^rcm-^ 

will have cooling time shorter than the dynamical time and will soon cool to T ^ 10^ °K 
, thus satisfying the Jeans instability criterion. The two conditions (4) are necessary and 
sufficient conditions for selecting gas particles as prone to SF. For the local SFR, NW adopted 
equation (2) with Cy^ = 1, = Td and £i, = 1/2 as the condition for which a gas particle 
can convert part of its mass into a star particle. These two algorithms will be referred to 
hereafter as KWH and NW, respectively. 

The numerical tests clOO — k and cl39 — k have been performed following the NW 
prescription for selecting gas particles which can form stars. The NW method has been 
preferred over KWH because of its simpler assumptions about the physical conditions of the 
gas in star forming regions. Because of the many physical processes involved in SF, having 
a minimal number of assumptions can reduce possible biases in hydrodynamical simulations 
when modeUing the local SFR. For a single representative cluster, a detailed comparison 
has been made between the cluster properties obtained using the two methods and different 
input values for the SF parameters. These simulation runs with numerical modelling of SF 
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also include energy feedback to cluster gas from supernova (SN) explosions. Once a star 
particle is created it can release energy into the surrounding gas through SN explosions. 
This energy is converted into heat of the neighboring gas at each time step, according to 
the stellar lifetime and initial mass function. A standard Miller-Scalo (1979) mass function 
has been adopted in the mass range from 0.1 to IOOMq. All of the stars with masses above 
8Mq end as SNe, leaving a IAMq remnant. Each SN explosion produces Ssn — 10^^ erg 
, or ~ 7.5 • 10^^ erg /Mq, which is added to the thermal energy of the gas. Current time 
steps are much smaller than stellar hfetimes, and so the SN energy is released gradually 
into the gas according to the hfetime of stars of different masses. At each time step, the 
fraction of stars releasing their energy into the medium is calculated for any star particle and 
the corresponding SN energy is spread over neighboring gas particles according to the SPH 
smoothing prescription. SN explosions also inject enriched material into the intracluster 
medium, thus increasing its metallicity with time. According to Steinmetz & Muller (1994) 
pz = 0.357m — 2.2 solar masses of heavy elements are synthesized by a SN progenitor of mass 
m. The enrichment in metals of the intracluster medium is modelled as follows (Steinmetz 
& Mueller 1994; Carraro, Lia & Chiosi 1998). Each SPH gas particle initially has zero 
metallicity; star particles are produced with the metallicity of the parent gas particle at the 
epoch of their creation. The metallicity of gas particles is successively enriched at each time 
step according to the fraction of exploding SNe associated with each star particle. The mass 
in metals produced by these explosions is calculated in accordance with the specified function 
Pz, and is added to the metallicities of the gas neighbors of the star particle. This mass is 
distributed over the neighbors using a smoothing procedure identical to that implemented 
for spreading the SN feedback energy of star particles among internal energies of the gas 
neighbors. According to the same procedure, the current mass fraction of exploding SNe 
with M > SMq is also added to the mass of the gas particles, with the exception of the 
1.4M(ri remnant. 



3. Results 

3.1. Cluster simulations with radiative cooling and cooling plus star formation 

The radial density and temperature profiles for the pure cooling runs are shown in Fig- 
ures 1 & 2. The cluster center has been identified as the maximum of the gas density. For 
each radial bin spherical averaged quantities have been obtained by estimating hydrody- 
namical variables at 100 grid points uniformly spaced in angular coordinates. Densities and 
temperatures at the grid points were computed from SPH variables according to the SPH 
smoothing procedure. The cluster ACDMOO is a particularly neat example of the effects at 
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work in the simulations. As can be inferred from the softening values reported in Table 2, the 

numerical strategy has been first to run clOO — 00 with the additional cooling prescriptions 
(clOO — 01) and in subsequent runs the value of Eg has been reduced in order to resolve the 
core radius of the gas density profile. Figure la shows that this is not achieved: whatever 
is the value of Eg, there is no evidence of a gas core radius, the gas density continues to rise 
steeply at the cluster center without any indication of converging to a constant value. This 
result is in strict agreement with those of others (Suginohara & Ostriker 1998; Anninos & 
Norman 1996; Pearce et al. 2000; Lewis et al. 2000) and it is known as the overcoohng 
instability; at the cluster center cooling times are very short because the gas density is high, 
thus drawing in more more material. The fact that the gas central density continues to 
rise as the spatial resolution is increased suggests that the extent of the (physical) effect is 
limited by the numerical resolution of the simulation. Of the four test clusters, CHDMOO 
is the only one which does not show the cooling instability with the exception of clOO — 05 
(Figure 2a). The reason for this behavior may be dynamical. Buote & Tsai (1996) measured 
a negative correlation of the cooling rate with the cluster X-ray substructure. The cooling 
instability can then be strongly suppressed when the cluster is still in a young dynamical 
state, with a large fraction of substructure. For the cosmology considered, VGB found that 
the clusters studied in the CHDM model with fi^ = 1 had many more substructures than 
those in the low-density ACDM model. In this case CHDMOO could be marginally stable, 
with instability being triggered by numerical effects when the central value of the gas density 
increases because the simulation resolution is increased. 

In the simulation runs the number of particles was increased as Eg was reduced so as to 
keep 2-body heating under control. The relaxation time , due to 2-body effects, is defined 
as 



where ai is the 1-D dark matter velocity dispersion, G is the gravitational constant, pd 
is the dark matter density; In A is the Coulomb logarithm, with A ~ Rh/AEg, and Rh is 
the half-mass radius. Typical values are In A ~ 3; standard theory gives A ~ Rh/^-, the 
factor 4 above accounts for the softening bimodal distribution (Farouki & Salpctcr 1982). 
For the simulations clOO — 00 and cl39 — 00 the relaxation time has been estimated at 
radius ~ 0.05r2oo — lOOi^pc, approximately the resolution limit of these simulations. At this 
length scale pd/ Pc — 2 • 10^. The values of range from ~ ISGyr (clOO — 00) to ~ VJGyr 
(cl39 — 00). For the simulation runs with increased resolution the minimum resolvable scale 
is set by Eg, and p/pc — 5 ■ 10^ at a fiducial scale ~ SOXpc. These rather high values of pd 
are a consequence of the slightly steeper profiles for dark matter in the inner regions, with 
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~ 6.7 • IQ^Gyr ( 




X 



(md/lQiiMo) (pd/pc)lnA ' 



(5) 



- 12 - 



respect to the case with no coohng. For these runs, the increase in is compensated by a 
corresponding increase in the particle number and hence a smaller value for m^, so that 
is approximately constant at the scale considered and is close to the Hubble time. Another 
timescale which is relevant for these simulations is the cooling time Tc, defined as 



where is the Boltzmann constant, n is the gas number density, T is the gas temperature 
and Ac ~ 5.2 ■ 10~^^T^/^n^erg'sec~^cm"^ is the cooling function. In the central gas regions 
Tc « Hq^ and a cooling instability will develop. In Figure 5 Tc is shown as a function of 
radius for the four test clusters. For the coohng runs clOO — 05, cl39 — 05 Tc is always below 
~ 20Gyr for r ^ lOOKpc, with the exception of ACDM39 , which has a bump but then a 
strong fall in Tg proceeding inwards. For the sake of reference, Tc for the no-cooling case 
( —00) has also been plotted. Because of the presence of a gas core radius in this case Tc 
approaches a constant value towards the center. According to Steinmetz & White (1997) 
gas cooling will be affected by artificial 2-body heating unless Tc(r) < Tr{r). This condition 
is satisfied if the dark particle mass is smaller than the critical value 

Me = 2 ■ lO^Te/cosMo, (7) 

where Tq is the the gas temperature in units of 10^ °K, /0.05 is the ratio / = Pg/pd in units 
of 0.05. For the simulated clusters studied here Te ~ 50 — 100, /0.05 ^ 0.5 for r ^ lOOKpc 
and Mc is always above ma. The simulations can then be considered free from numerical 
effects which can dominate the gas behavior. In the cooling simulations this condition might 
be violated near the cluster center, where Tq ^ 10 for r ^ 50Kpc. However, in these regions 
Tc << Td — 27h~^Gyr / ^J~Pdfpc and the cooling is effective in removing the gas energy at a 
faster rate than the one set by dynamical effects. 

The temperature profiles show a decrease for r ^ lOOXpc and a drastic drop in the 
central values, where cooling is most effective. This inversion in temperature takes place in 
all of the tests considered. Peak values for the gas temperature are located at ~ lOO/Tpc 
(0.05—0.1 of r2oo)- Between this distance from the cluster center and r2oo the gas temperature 
decreases with radius and the clusters clearly cannot be considered isothermal. These results 
are in agreement with those of Pearce et al. (2000), and suggest that the global cluster 
properties are affected by coohng processes active on inner scales , where the coohng time is 
short (see also Lewis et al. 2000). The most important cluster variable which is affected by 
these results is the cluster X-ray luminosity. For evaluating Lx the standard SPH estimator 
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gives 

1 

Lx = 5.2 • 10-28.^-—^ J2 miPiTl^'ergsec-\ (8) 

where, for the coohng function, the bremsstrahlung emissivity has been approximated with 
a gaunt factor of 1.2 (eq.[2] of Suginohara & Ostriker 1998), // = 0.6 is the mean molecular 
weight, and nip is the proton mass. The summation is over all of the gas particles within 
r2oo- Figure 7 shows the behavior of Lx at z = as a function of Eg for the pure cooling 
runs (open symbols). There is not a clear convergence of Lx as the resolution is increased. 
In fact Lx obeys the approximate scaling Lx oc l/^^^- Similar results for Lx have been 
obtained by Anninos & Norman (1996) in their convergence study of simulations of X-ray 
clusters. 

The unphysically high values found for Lx in the pure cooling runs arise because the gas 
density continues to increase steadily at the cluster center, while the conditions of high gas 
density and low temperature cause the gas to become Jeans unstable. Thus the treatment 
of gas cooling in cluster simulations cannot be decoupled from a modelling of the physical 
processes turning the cold, dense, gas into stars. The coohng simulations have therefore been 
rerun with the inclusion in the integrations of an algorithm for converting gas into stars. 
This used the NW method, with parameters = 1 and £^ = 1/2. In these simulations 
gas particles can produce star particles without any limit on the number of star-forming 
events. Furthermore, the SN explosion energies and metallicities that star particles can 
produce are smoothed over 32 gas neighbors but with an upper limit of Hm = 15Kpc for 
the SPH smoothing length. This is in order to avoid unphysical heating of the gas over 
length scales much larger than those involved in the SF activities. This upper limit is also 
justified by the lack of diffusion in the ICM of the metals injected from galaxies (Ezawa et 
al. 1997). The results obtained are shown in Figures 3 and 4; the index of the simulations 
is A; = J -|- 5, where j is the index of the cooling runs. The most important result is that 
the inclusion of an SF model has been effective in removing the unphysical gas behavior, 
which now shows a well defined core radius in the radial density plots. This is valid in all of 
the cases considered, with the expection of CHDMOO. The simulation runs clOO — k do not 
show, for this cluster, evidence of an SF activity, a result which is in agreement with the lack 
of a cooling instability in the cooling simulations clOO — j. The shapes of the temperature 
profiles show that convergence is achieved for Ng ^ 20, 000. For the simulation runs with 
the highest resolution, all of the central values for the gas temperatures at r = 10 Kpc 
are within a factor ~ 1.5. CHDM39 is an exception to this rule, with cl39 — 10 still not 
showing a flat temperature proflle in the inner regions and resembling that of cl39 — 05. It 
is unlikely that the source of the discrepancy is due to a convergence problem: cl39 — 11 has 
a slightly larger accuracy in the computation of the gravitational forces (Hernquist 1987), 
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nevertheless its temperature profile has these discrepancies largely removed. For CHDM39 
the peculiarity of cl39 — 10 in the final temperature profile, with respect to the other numerical 
tests shown in Figure 4b, could be of a numerical nature: for a certain accuracy in the tree 
evaluation of the gravitational forces, matter subclumps might form during the integrations 
which can then modify the gas dynamics. The formation of these sub-clumps is triggered 
by the approximations involved in the truncation of the multipole expansion of the cluster 
gravitational potential. The statistical occurrence of this effect should be small, because it 
is not observed in the other three test clusters. For a tree method the errors involved in the 
multipole expansion of the gravitational potential have been estimated by Hernquist (1987) 
assuming a spherical, isotropic Plummer model for the mass distribution of N test particles. 
A comparison against the accelerations obtained by a direct sum shows that in the large 
N limit (N ^ 30, 000) the errors in the tree evaluation of the accelerations are negligible 
for a monopole expansion if ^ 1. The inclusion of the quadrupole terms improves the 
accuracy of the force computation, for 9 1 the errors in the forces arc those of a monopole 
expansion with 9 0.8. The convergence in T(r) obtained for cl39 — 11 suggests that for a 
given accuracy, quadrupole corrections should be preferred when evaluating tree forces. 

The radial density profiles of the star component are also shown for the various runs in 
the density plots. The slope of these profiles is approximately ~ —3, a value close to the one 
observed for galaxy populations in galaxy clusters. In all of the simulations, the gas density 
profiles have a well-defined core radius, with size Tc — 50 — lOOKpc, approximately 0.05 of 
the virial radii. From the density profiles note also that the gas core radii are smaller than 
in the no-cooling runs, outwards of Tc the density profiles are very similar to the no-cooling 
cases. The temperature profiles increase inwards from the virial radius up to ^ l(}(} — 200Kpc. 
Thereafter the profiles stay almost flat, or with a modest decrease in T(r) towards r = 0. 
The strong drop of the temperatures in the very central regions for the cooling runs is no 
longer seen, the inclusion of a star formation prescription having been effective in removing 
the cold gas particles ( ^ 10^ °K) from the chister centers. 

Cooling timescales Tc(r) are plotted in Figure 5; the dashed line in the four panels is 
for the simulation runs including SF. As a general rule, for each test cluster, the Tc are 
almost indistinguishable in all of the simulation tests for r ^ lOOi^^pc. In the simulations 
including SF, Tc moves toward the no-cooling case in the cluster inner regions because of 
the reduction in the gas central density. The central values of Tc are well below the present 
age of the universe for all of the models, with the exception of CHDMOO. This cluster does 
not show a coohng instability and has Tc ~ 20Gyr in the cluster core. Accretion rates 
M(r) = 4:7rpgr'^Vr are plotted in Figure 6 for the same test clusters as in Figure 5. For each 
radial bin, spherical averages of M(r) are shown only for negative values of f ,.. For the test 
clusters in the A-dominated cosmology there is a well defined radial infall of matter within 
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r ^ lOOKpc, compared to the adiabatic run. These inflows of matter can be compared with 
those estimated from X-ray data for coohng flow clusters. For example Thomas, Fabian 
& Nulsen (1987) present mass deposition profiles M(< r) for a sample of 11 cooling flow 
clusters. They obtain for Abell 478 M ~ 10^ MQyr~^ at r ~ SOOKpc (Fig. 7 of their paper). 
For this cluster they quote a measured line-of-sight velocity dispersion of ~ 750Kmsec~^ 
and the estimated cluster virial mass ( ~ 4 • lO^'^h'^MQ,) is close to that of ACDM39. This 
accretion rate is in good agreement with the values shown in Figure 6 at r ~ SOOKpc for 
ACDM39 in the simulation including SF. Note that in order to compare M{< r) with M(r) 
one is implicitly assuming a steady-state. The clusters in the cosmology with = 1 show 
values of M which are higher in the adiabatic case, in comparison with the simulation runs 
including cooling and star formation. The reason for this discrepancy is of a dynamical 
nature: CHDMOO, for example, has a large radial infall velocity because it is still out of 
equilibrium, with material collapsing onto the center. These results seem to support the 
hypothesis of an anti-correlation between the strength of the cooling flow and the dynamical 
state of the cluster, as measured by the power ratios, for example. Buote & Tsai (1996) 
demonstrated the existence of such an anti-correlation, as expected from a dependence of 
the cooling flow rate on the cluster dynamical state. A statistical analysis of this correlation 
is beyond the scope of this paper and is left to future work, where the analysis will be 
performed for the whole sample of 40 clusters for each of the three cosmological models 
(VGB). A striking result is the convergence of the final X-ray luminosities for the simulation 
runs including SF for the four test clusters. In Figure 7 the filled symbols refer to these 
simulations. The divergence for Eg ^ of the cooling runs is completely removed and the 
plotted values are quite stable. An exception is cl39 — 10 of CHDM39 (the black square 
in the bottom right panel) which has Lx — 5 • lQi^ergsec~^ . As previously discussed, the 
peculiarity of this cluster is of a dynamical nature and there is not a question of convergence 
of hydrodynamical variables. In fact cl39 — 11 has Lx — 2 ■ 10^'^ ergsec~^, a value in full 
agreement with the values of the other runs. Compared with the non-radiative case, the 
luminosities are stable (CHDM) or increase by a factor ~ 2 (ACDM). 



3.2. Comparison with previous simulations 

These results can be compared with previous findings from hydrodynamical cluster 
simulations, in which the gas is allowed to undergo cooling and star formation or is subject 
to a prescription for the treatment of cold particles. The density plots can be compared 
to analogous plots shown in Lewis ct al. (2000). In their paper Lewis et al. analyzed five 
simulations of a cluster with M200 — 4- 10^"^ M© in a CDM universe with Qm = 1 and h = 0.5. 
Of these simulations, the Cool + SF allows gas to undergo cooling and star formation, using 
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the KWH method with = 0.1 and Esn — 1- Although the cosmologies are different, a 
rough comparison can be made for ACDM39, which has the closest virial mass to their test 
cluster. For ACDM39 the simulations cl39 — 10 or cl39 — 11 have a numerical resolution 
comparable to the Lewis et al. (2000) simulations (compare Tables 3 & 4 with Table 2 of 
Lewis et al.). Thus the density plots in Figure 3b can be compared with Figure 7 of Lewis 
et al. . The results are encouraging, there is a rough agreement for the various baryonic 
components, although in Lewis et al. there is a central spike in the gas density which is 
not observed in the present simulations. A substantial difference is instead found for the 
temperature profiles: all of the simulation runs show a tendency for T(r) to recover an almost 
fiat behavior in the inner regions (with the exception of cl39 — 10 for CHDM39, which is 
peculiar for the reasons previously outlined). Lewis ct al. found instead that T{r) in their 
Cool + SF simulation reaches a peak value of ^ 8 • 10^ °K within ^ 40Kpc (Figure 9 of Lewis 
et al.). The radial behaviors of Figure 5 (Tc(r)) compare well with Figure 8 of Lewis et al. 
(2000). 

Hydrodynamical simulations of cooling clusters have also been analyzed by Pearcc ct 
al. (2000) and Yoshikawa, Jing & Suto (2000). I discuss here in detail how the results of 
§3.1 compare with the cluster properties of Yoshikawa, Jing & Suto (2000, hereafter YJS). 
The simulated clusters of the Pearce et al. (2000) runs have properties similar to the YJS 
clusters, the only substantial differences being found for the X-ray luminosities. 

In their paper YJS analyze results from a set of cosmological SPH simulations and 
concluded that estimates of the X-ray luminosities arc biased by the numerical resolution of 
the simulations and are not reliable. The cosmological simulations of YJS were performed in 
a flat, A-dominated CDM cosmology. The values of the background cosmological parameter 
are identical to the ones chosen here for the cluster simulations in the ACDM model. They 
performed simulations with two different box sizes : L — 75h~^Mpc and L — 150h~^Mpc. 
They took 128^ gas particles, an equal number of dark particles, and a comoving softening 
parameter e = L/1280. In the simulations including radiative cooling, cold gas particles were 
removed from the summations defining local gas variables if the Jeans condition in equation 
(1) is satisfied (simulation runs with label UJ in their notation). This constraint is slightly 
different from the one of Pearce et al. (2000), but is almost identical for Tj ^ 10"^ °K. Their 
prescription is phenomenological and is intended to take into account the process of galaxy 
formation. For the cosmological simulation with box size L — 150h~^Mpc, the numerical 
resolution for the most luminous clusters is comparable to that of the clOO — 11 run. The 
most massive cluster in the simulation with label 150UJ has ~ 1.56 ■ IO^^Mq, a value close 
to that of M200 in Table 1 for ACDMOO. For this cluster YJS found at z = a bolometric 
X-ray luminosity ~ 2 • W^^ergsec'^ , about twice the corresponding value of Lx for clOO — 11. 
A comparison between the cluster density and temperature profiles (Figure 1 of their paper. 
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central panel) and the analogous plots of the present paper (clOO — 11 in Figure 3) shows 
that there are substantial differences at r ^ bOKpc. In YJS the cluster gas density shows 
no evidence for a core radius; the mass-weighted temperature profile has an inversion at the 
cluster center, as also found by Pearce et al. (2000), while clOO — 11 shows a nearly fiat 
profile for T(r). These differences in the gas profiles at the cluster centers are probably due 
to the different methods employed for the treatment of the cooled gas, rather than to the 
numerical resolution of the simulations. 

YJS adopted the phcnomcnological prescription of treating separately those gas particles 
which satisfy the Jeans criterion. Therefore the gas distribution consists of gas particles in 
a 'hot' X-ray emitting phase together with a population of cold particles at temperatures 
around 10^ °K, the latter being locahzed at the cluster center. The temperature profiles in 
YJS include the cold gas population ^ and show a steep decrease at the cluster center. In 
the simulation run clOO — 11 this is not observed because the star formation prescription 
has been effective in converting the cooled gas at the cluster center into the form of stars. 
For the simulation runs with cooling and star formation, the mass-weighted temperatures 
Tm{sim) of the clusters ACDMOO and ACDM39 can be compared with those expected from 
the isothermal mass-temperature relation. At the present epoch, this relation reads: 



1/3 / , . N 2/3 



ksTx ^ 5.27 —-f- KeV, (9) 

where M is the cluster virial mass. A,, — 178Q^'^^ is the virializcd cluster ovcrdensity and 
7 is a fudge factor. YJS found that the mass-weighted temperatures of their simulated 
clusters are weU fitted by equation 9 with 7 = 1.2. For the cluster ACDMOO (ACDM39) 
I find from the simulations clOO - ll(cl39 - 11) that T^{sim) = 5.6{3.1)KeV, while the 
theoretical relation (9) gives T^{th) — b.^{2.8)KeV . Thus the mass-weighted temperatures 
of the simulated clusters are in close agreement with the theoretical predictions and also 
with those found by YJS in their simulations. 

As previously shown in the simulation runs including star formation, the X-ray luminosi- 
ties are found to be numerically stable and converging to reliable values. With respect to the 
adiabatic runs, the Lx are found to be constant or with a modest increase (Figure 7). These 
results are at variance with those of Pearce et al. (2000), who found that cooling clusters are 
less luminous than those in the no-cooling runs. Similar results have been obtained by YJS 
and Lewis et al. (2000), who measured an increase in the final X-ray luminosity when cooling 
was included. A comparison with the Cool-I-SF simulation of Lewis et al. (2000) is difficult 
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because of the different cosmologies and cluster virial masses. YJS analyzed the Lx — Tx re- 
lation obtained from the simulated clusters in the UJ test runs at ^ = 0, from the simulations 
with two different box sizes. The cluster luminosities in the simulation with L = 150h~^ Mpc 
are found to be underestimated with respect to the ones of the L = 75h^^Mpc simulation 
box. YJS draw the conclusion that, in order to have reliable cluster luminosities, simula- 
tions of cooling clusters require a much higher numerical resolution than the one employed 
in their simulations. In the simulations with L — 150h~^Mpc, the mass of the gas particles 
is ~ 2 • W^Mq, a factor ~ 2 larger than that of clOO - 10 for ACDMOO. Therefore their 
conclusion seems to be in conflict with what is claimed here, that substantial convergence 
in X-ray luminosities is achieved for the highest resolution simulation runs. The source of 
this discrepancy lies in the different numerical approaches. YJS simulated cluster evolution 
in a cosmological box with a constant gas particle mass. As outlined by the authors, the 
resolution problem is severe for the less luminous clusters in the simulation box. The multi- 
mass technique described in §2 is instead used here to simulate a single cluster. The results 
of the numerical simulations show that convergence in the gas variables is obtained for each 
single cluster whenever Ng ^ 20, 000. A comparison with the numerical resolution adopted 
by YJS is useful. In their simulations with L — IbOh'^Mpc the gas particle mass is com- 
parable to the value found here for which ACDMOO can be safely analyzed. However they 
have a constant mass resolution and their cluster sample has a lower limit of M > IO^^Mq. 
An application of the numerical parameters required here to a cluster with a virial mass of 
~ 10^'' Mq would give a gas particle mass of ~ ^li.lO^'^ / {ilm2.2 ■ 10^) ~ 5 • 10^ Mq, a factor 
~ 30 smaller than that being used in the L = 150h~^Mpc simulations. These values also 
show that cosmological simulations require a number of gas particles ~ {AOOL / 150h~^ MpcY 
in order to give realistic estimates for the statistical properties of X-ray clusters. 



3.3. Simulations with different star formation prescriptions 

The numerical tests studied here give the range of numerical parameters for which the 
results of cooling cluster simulations, including SF according to the NW prescription, reach 
numerical convergence and can be considered stable. A different question concerning the 
reliability of the numerical results is the sensitivity of the estimated cluster properties to 
the numerical method used to describe star formation and energy feedback from stars in the 
hydrodynamical simulations. To investigate this, final results for different simulations have 
been compared for a single test cluster (ACDMOO). The simulations were performed with 
the same numerical parameters as for the clOO — 11 run, but using different SF methods or 
parameters. This was done in order to demonstrate the effects of the SF modelling on the gas 
variables. A summary of these simulations with different SF prescriptions is reported in Table 
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6. The simulation with index I is the standard case with which previous coohng+SF runs 
have been performed. Thus this simulation just corresponds to cIOO — 11. Two simulation 
runs correspond to the NW method but with a different feedback energy for SN explosions 
(1 and V). In the other three runs the KWH prescription was adopted for converting gas 
particles into stars. Two of them (II and III) compare the results obtained for a different 
star formation efficiency parameter c*, with the other parameters being held fixed; in a 
third run (IV) the fraction £^ of mass of the gas converted to stars is varied. The results 
of the different approaches are shown in Figures 8 & 9. For the simulation V, only the star 
formation rate (SFR) and the X-ray luminosity versus time have been plotted in the two left 
panels of Figure 8. This is because the simulation V produced final results almost identical 
to the reference case I. In simulation V the feedback SN energy was set to Esn — lO^'^erg, a 
value ten times smaller than that of the simulation run I. A comparison between the results 
plotted in Figure 8 shows that the two simulations have an identical evolution for the X-ray 
luminosity, but a different SFR. This is because of the smaller amount of SN feedback energy 
which is added to the thermal energy of the gas in the run V. This in turn implies higher 
gas densities and SF rates for simulation V. The differences become negligible for t ^ SGyr. 
The final gas profiles, as well as the other variables, are identical. 

The simulations with the KWH method and different ( II and IV), give similar 
results and show that the choice of the mass fraction e^, is not important in modelling the 
star formation processes. The most important differences are found between the KWH 
simulations with different (II and IV). The differences are dramatic in the final X-ray 
luminosities, which differ by a factor ~ 40. The source of this discrepancy hes in the 
different gas density profiles, which have substantial differences in the cluster core regions 
for r ^ lOOKpc. These differences are localized at the cluster center; beyond r ~ lOOKpc all 
of the profiles converge, as shown in the plots of Figure 8. The temperature profiles have a 
peak value of ^ 10^ °A' at r ~ lOOKpc and thereafter decline outwards by a factor ~ 2 out to 
r2oo- Below ~ lOOKpc the profiles instead show large differences. Compared to run I (NW) 
the two simulations with = 0.1 have gas temperatures which decrease inwards by a factor 
~ 10 from ~ lOOKpc down to r ~ lOKpc. These radial decays follow because of the less 
efficient conversion of the cooled gas into stars compared to the NW run. These results can 
be compared with those of Lewis et al. (2000), who used the same star formation method and 
parameters, although in a different cosmology. The temperature profile of their Cool+SF 
simulation shows a similar decrease outside of the cluster core (Fig. 9 of their paper). Is is 
difficult to compare the radial runs at r < lOOKpc because the authors adopt a linear scale 
for the plots. The gas density profile of run II rises steeply towards the cluster center from 
r ~ lOOKpc; a similar behavior is shown by Lewis et al. (2000) for their hot gas population 
( see their Fig. 7 ). The simulation run with Ci, — 1 has radial profiles much closer (but not 
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identical) to the NW ones. These differences are also reflected in the estimated emission- 
weighted temperatures (see Table 6). For simulations I and III — 8 — 9KeV, while 
Tern — 2KeV for run II. Note that for this simulation most of the contribution to T^,^ comes 
from within the 50Kpc cluster core. 

There is a remarkable agreement for the ratio of the cluster mass locked into stars to the 
gas mass, which is ~ 10% at r2oo for all of the runs considered (Figure 9, bottom left panel), 
with respect to the observational values ( Evrard 1997). Also the density profiles of the stars 
produced are very similar. As a general rule, one can say that the global cluster properties 
are not affected by the SF prescription adopted. The soiuce of the higher gas densities in 
the cluster core for the simulations with = 0.1, compared to what is found for = 1 , 
is in the different SFR that the two simulations use during the integration. In the top left 
panel of Figure 8 the SF rates are plotted for the different runs. The rates are plotted versus 
time instead of redshift because otherwise final differences would have been compressed. It 
can be seen that the simulations II and IV have a different SF histories, with respect to runs 
III and I. At early times runs II and IV have a lower SFR than the simulation with = 1. 
This behavior is reversed after t ~ lOGyr [z ~ 0.3), with simulations II and IV showing a 
substantial SF activity, while runs III and I have a sharp decline in their SF rates. As final 
results, the two simulations II and IV have less gas converted to stars in the cluster core. 

According to KWH, in a simulation with = 0.1, the collapsing gas will reach an 
equilibrium density higher than in the one with c+ = 1. This is because of the reduced 
efficiency in the conversion of the gas into stars for the simulation with c+ = 0.1, which 
implies that gas collapse will proceed until pressure from the thermal energy of the gas is 
able to prevent further gas collapse, reaching an equilibrium at higher gas densities. This 
in turn implies higher SF rates so that, after an initial transient, simulations with c^i, = 0.1 
will have higher gas densities and SF rates compared to the — 1 runs. This expected 
behavior is what is found in the plots of Figure 8, the NW simulation giving identical results 
to the KWH run with = 1. The bottom right plot of Figure 9 shows Ms{z) the cluster 
mass in stars at the redshift z and is particularly useful in analyzing differences in the SF 
histories between different runs. The mass in stars for the NW run stops growing aX z = 0.5, 
and stays flat until z = 0. The KWH simulation with = 1 has a similar trend but 
with a somewhat steeper slope in the growth of Ms{z). An important result is that the 
final values of the cluster mass in stars vary by only ~ 20% from Mg ~ 1.5 ■ IO^^Mq. As 
previously outlined, this shows the reliability of different SF prescriptions in reproducing 
global cluster properties. In the simulation results, differences between the methods arise 
in the SF rates, which are determined by different choices of the star formation efficiency 
parameter c*. In Figure 9 for simulations II and IV, the values of Ms{z) at high redshifts 
are below the corresponding values of the c* = 1 runs. For these simulations with c* — 0.1, 
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Ms{z) has a continuous growth and is smaller than Ms{z) of the NW run until z ~ 0.2. 
Note that the final masses in stars for the different runs do not correspond to what would 
naively be expected from the different shapes of the gas density profiles in the cluster core, 
Ms{z = 0) being determined by the overall SF history. The SF rates of the simulations 
with Cv, = 0.1 have the expected behavior, with respect to those of the c^, = 1 runs, but the 
high resolution runs of the next section will show that these SF rates are depressed by the 
numerical resolution of the simulation. 

The emission weighted metallicity profiles at z = are shown in Figure 10. Conversion 
from the mass fraction in metals associated with each particle to the iron mass was accom- 
plished according to the relationship between the total metallicity and the iron abundance 
[Fe/H] of Tantalo, Chiosi & Bressan (1998). The relative iron abundances Z — Fe/H in 
the radial bins were estimated from the individual iron mass of each gas particle using the 
SPH smoothing procedure. These profiles were still rather noisy, and to obtain the final 
radial profiles of iron abundance, a further smoothing was performed by considering only 
five distinct radial bins and averaging over neighboring bins. An important result for the 
metallicity profiles is that the simulation runs show the existence of radial gradients, with 
decreasing metallicities as the radius increases. This is in broad agreement with observa- 
tional data for cooling fiow clusters (Ezawa et al. 1997; Kikuci et al. 1999; Buote 2000; 
White 2000; De Grandi & Molendi 2001). For the simulation runs I and III Z > 0.5Zq at 
the cluster center and Z ^O.IZq at r ^ lOOKpc , with a value for the solar abundance 
Zq = {Fe/H)Q = 4.68 • 10"^ (Anders & Grevesse 1989). At r > lOOKpc these values 
are ~ 1/3 of the measured abundances obtained for a sample of 9 cooling fiow clusters (De 
Grandi & Molendi 2001). This deficit of iron abundance is probably due to the lack of metal 
enrichment of the intraclustcr medium from SN of type la and will be analyzed in a future 
paper, where the metallicity dependence of the cooling function will be implemented in the 
simulations. The KWH runs with = 0.1 have metallicity profiles with abundances well 
below the NW run, this is because of the different SF histories, with simulations 11 and IV 
having a lower SFR at early times ( ^ IQGyr). The KWH run with 0^,-1 has, on average, 
a steeper metallicity gradient than the NW run; here again this is because of the different 
SF histories, as shown by the different growth of Ms{z). 



3.3.1. High resolution simulations with different star formation prescriptions 

An important point about the simulation results related to the SF parameters, such 
as for example the SFR, is their dependence on the numerical resolution adopted. The 
numerical tests of Table 6 have been performed with the same number of particles and 



-22- 



softening parameters (clOO — 11 in Table 3), and for the particle gas mass rUg ~ 1.2 • lO^^M©. 
The resolution of the mass of gas is clearly inadequate for modelling a single galaxy formation 
process. SPH simulations of galaxy formation processes have been debated by many authors 
(Thacker & Couchman 2000, and references cited therein). In order to resolve the internal 
dynamics and to follow shock evolution of a forming galaxy a minimum of ~ 10"^ particles 
is required (Thacker et al. 2000). Lia, Carraro & Salucci (2000) discussed the dependence 
of the mass resolution on the SFR in gas-dynamical simulations of a collapsing spheroid. 
According to their results, the total mass of gas converted into stars diminishes from 90% to 
~ 80% of the initial mass of gas, passing from the low resolution run with 2, 000 particles to 
the high resolution run with more than ~ 10^ particles. However, early simulations (Evrard 
1988, NW) showed that SPH simulations with even a small number of particles (~ 100) can 
converge to stable results as far as global properties are concerned. For example Figure 20 
of NW indicates that the final mass of stars in their SPH simulations of a rotating cloud is 
already close to the convergence value for a number of gas particles Ng ^ 100. 

It is therefore important to assess the effects of numerical resolution on final results 
in the simulations with different SF prescriptions. This has been done by running again 
simulations I, II and HI of Table 6 but with a number of particles increased by a factor 
~ 3. For the gas particles Ng = 69,599, rUg = 2.4 • lO^/i^^M© and Eg = lO.bh-^Kpc. The 
cold dark matter particles have these values scaled in proportion. These simulations will be 
referred as IH, IIH and IIIH, respectively. The numerical parameters are those of clOO — IIH 
in Table 3. Different SF parameters correspond to the ones of Table 6. The simulation 
results are shown in Figure 11. For simulations IH there are not appreciable differences 
in the radial profiles. This confirms the results of §3.1, that the NW runs have reached 
numerical convergence in the physical variables for the numerical parameters of the clOO — 10 
simulation of Table 3. The profiles of simulation IIH are instead different from those of run 
II at r ^ 50Kpc. The strong drop in T(r) has been removed and the gas density profile is 
much closer to the NW one. Simulation IIIH gives final profiles very similar to the ones of 
the parent simulation. The bottom left panel of Figure 11 shows that high resolution runs 
have final X-ray luminosities which can differ by a factor ~ 2 from the parent simulations. 

SF rates are shown in the the top left panel of Figure 11 and at 2; > there arc large 
differences between the high resolution run IIH and the parent simulation. For t ^ 5Gyr, 
simulations IH and IIIH give the best performances, with the SF rates closely following the 
ones of the lower resolution runs. Note the strong dechne in the SFR after 2; = 0.3 . At 
early times {t ^ 5Gyr) the SF rates of simulations I and HI are much lower than those of 
the corresponding high resolution runs. For these runs, the peak of the SF activity is shifted 
from z ~ 0.7 up to ^ ~ 2 (t ~ 3Gyr). This shows that in order to correctly sample the 
SF rate over the whole cluster evolution the numerical resolution must be at least that of 
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the high resolution runs. For t ^ 5Gyr, simulations I and III have cluster SF rates in good 
agreement with those of the corresponding high resolution runs, while simulation II does not 
show that convergence is achieved when the resolution is increased. For simulations I and 
III, the X-ray variables are not affected by the undersampling of the SF rates at early times; 
for example, the X-ray luminosities (Figure 11, bottom left panel) are fairly stable under an 
increase in the numerical resolution. 

The results of these high resolution simulations show that in the large A''^ regime, differ- 
ent star formation methods approach similar gas profiles. For the simulations of Table 6, the 
differences found between the profiles of the runs with = 0.1 and the others arise because 
in the former case, with respect to the KWH run with = 1, the increase in the SF rate 
corresponding to a decrease in the value of is limited by the numerical resolution. The 
upper left panel of Figure 11 shows that for the high resolution run IIH the SF rate is much 
higher than in simulation II. 

From the simulation results it follows that the NW method is the most efficient in the 
removal of cold gas from the cluster center. As already stressed, in principle all the methods 
converge to the same profiles when Ngas gets very large. The runs with = 0.1 are the ones 
with the most important differences in the temperature profiles, compared to the other runs, 

and the reasons of these differences have been previously discussed. An explanation of why 
the NW method is so efficient relics on the chosen criterion for selecting cold gas subject to 
star formation, the method being based on a density threshold. If the local gas density pi 
exceeds this threshold (7.10~^^(?rcm^"^) then a mass fraction £^ = 1/2 of the gas particle is 
converted into a star particle. In very high density regions, the timescales of star formation 
are very short and the SF process removes the gas very quickly. Because the criterion is 
based on a density threshold this means that all the gas particles above this threshold are 
selected for a star forming event. 

Differences in the temperature profiles between the NW and KWH runs with = 1 are 
minimal and are localized at the cluster center, these differences being due to the different 
dependence of the two methods on the resolution limit of the simulations. The two criteria of 
SF depend on the numerical resolution of the simulations through the SPH smoothing lengths 
hi, which determine pi and are constrained by the lower limit hi ^ £g/4 = hmin, where Eg 
is the gravitational softening length. As the gas density increases toward the cluster center, 
the smoothing lengths hi get smaller, until they reach the limit hmin- The results of the 
simulations show that at the cluster center, in the small value regime hi — > hmin, the NW 
criterion for identifying cold gas particles is less sensitive to the resolution limit hmin than 
the KWH method based on the local Jeans instabihty. 

To summarize, the above results demonstrate that simulations I and HI of Table 6 have 
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an adequate numerical resolution to reliably predict X-ray cluster properties, such as the 
X-ray luminosity. For simulation II ( KWH with = 0.1 ) there are still differences at the 
cluster core between the final profiles when the numerical resolution is increased. Run IIH has 
a final Lx which is very large compared to the expected range of values from the luminosity- 
temperature relation (see below). For the simulation runs I and III Lx — IQ^^ er g sec~^ , 
while Lx — 4 ■ lO'^^ergsec^^ for the runs with c^, = 0.1. These large discrepancies can be 
reduced by adopting the phenomenological prescription of removing from the summation (8), 
cold gas particles with temperatures below a cut-off value Tc. For the runs II and IV with 
Ci, — 0.1, Lx — 10^^ ergsec~^ if the above prescription is adopted with Tc — 3 • 10^ °K. It is 
important to note that the NW run has an Lx that shows no sign of evolution for t ^ 5Gyr 
{z<1.2). 



3.4. Comparison with observational parameters 

These differences between different methods suggest that a rehable SF algorithm should 
be chosen by requiring that simulation results should consistently satisfy a wide set of differ- 
ent observational constraints on cluster properties. To this end, the results of the simulation 
runs reported in Table 6 have been compared with several cluster observations. The data 
investigated arc : the Lx — Tx relation for cooling flow clusters and the estimated SFR in 
rich clusters at the present epoch. These observations have been chosen because the plots 
of Figure 8 showed large differences between simulated clusters for the variables connected 
with these data. 

For cooling flow clusters, Allen & Fabian (1998a) have studied the relation between the 
bolometric X-ray luminosity L^x and the cluster temperature. They use ASCA spectra and 
ROSAT images to construct a sample of 30 luminous clusters {L^x > 10'^^ ergsec'^^), 21 of 
which have central coohng times < lOGyr and are identified as coofing fiow (CF) clusters. 
Their best-fit relationships are derived for a cosmology with — 1 and h — 0.5, thus the 
values reported in Table 1 of their paper must be rescaled to those for a fiat cosmology with 

— 0.3 and h — 0.7 . Allen Sz Fabian (1998) fit the L^x — Tx relation with a power-law 
of the form kTx = PL'^xj where kTx is the spectral fit temperature of a single isothermal 
model in KeV and L^x is in units of lO"^"^ ergsec^^. The data values used are those of model 
C in Table 1, which according to the authors gives the best values. Thus the coefficients 
P and Q obtained here for the A-cosmology arc those corresponding in Allen & Fabian 
(1998) to the only for the entry CFs in their Table 2. The least-square estimator gives 
P — 1.97 ± 0.33 and Q — 0.42 ± 0.05, error bars are the one-parameter 68% confidence 
limits, with x^/d-O-f — 1-03 and a mean gas fraction at 360Kpc of fg — 0.11 ± 0.03. For 
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the simulation runs I and III Lhx ~ 1.3 ■ W^^ergsec"^ and the predicted cluster temperature 
is T^'^'^ = 5.8 ± 1.2eV. The error bar represents the la confidence interval. The simulated 
cluster in the NW run has an emission-weighted temperature = 8.3KeV and a mass- 
weighted temperature = 5.6KeV. These values are consistent with those extracted from 
Figure 5 of YJS for the most massive clusters in their UJ simulations. They are also in 
good agreement with the the values predicted from the mass-temperature relation inferred 
from cluster ASCA and ROSAT data (Horner, Mushotzky & Scharf 1999; Nevalainen, 
Markevitch, & Forman 2000). For the simulation run III the values of L^x and the weighted 
temperatures are similar to the ones of the NW simulation. Mathiesen & Evrard (2001) 
have discussed how measurements of the intracluster medium spectral fit temperature Tg are 
related to and showed that Tg can be considered as a nearly unbiased estimator of T^. 
Thus the values of the mass-weighted temperatures reported in Table 6 can be compared 
with the corresponding T^'^'^. The simulations 1 and 111 are consistent at the \a level with 
the Lftx — Tx relation derived from a sample of cooling flow clusters. Simulations 11 and IV 
have Lfox — 4. ■ IQ'^'^ergsec^^ and the predicted cluster temperature is outside the 2a limits. 

Another observational test is to compare for the different simulation runs the rates of 
SF with the observed mean SFR of cluster galaxies. There is now a strong observational 
evidence that star formation in a cluster environment is strongly suppressed with respect 
that of field galaxies (Balogh et al. 1998; Poggianti et al. 1999). Kodama & Bower (2000) 
have estimated the SF histories in rich cluster cores for four clusters, using photometric 
models constructed from 7 CNOC clusters at 0.23 < z < 0.43, and found a strong dechne in 
the SF rate relative to the field for z < 1. The SF histories of the four clusters are shown in 
Figure 9 of Kodama & Bower (2000). One of these clusters is Coma, for which the measured 
values for the mass in galaxies, hot gas and cluster total mass (White et al. 1993) are in 
the same range as those of the simulated cluster ACDMOO. Thus a comparison between the 
SFR estimated from cluster simulations and the observed data can be made consistently. 
Although the numerical resolution of the simulations 1-111 is inadequate to correctly sample 
the SFR of the galaxies in the cluster, the previous discussion about simulations IH, IIH 
and IIIH has shown that even for the low- resolution runs of Table 6, the global cluster SFR 
is roughly reproduced. Therefore the SF rates computed in simulations I-III can be used 
for making a crude comparison with the values derived from Kodama & Bower (2000), in 
order to estimate the consistency of the simulation results for different SF algorithms with 
observations. 

The SF rates plotted in Figure 8 have been estimated by summing the mass of gas 
converted into stars in time bins of size Ai^ = 3 ■ 10^7/r. For simulations I and III, at 
z = Ats « Tg and the random process used to sample the distribution (3) can be 
approximated as Poissonian, thus the SFR ~ 70 ± 70h~'^MQyr~^; in fact there is only 
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one single event in the last four time bins. For simulation II, the condition Atg » Tg 
applies and the SFR ~ 330 ± 110h~'^MQyr~^; the dispersion has been computed over the 
last four bins. These values must be rescaled to the normalization adopted by Kodama & 
Bower (2000), who reported the integrated rates for galaxies in the cluster cores in units 
of 1O^^M0 within a radius enclosing 1/3 of the galaxy population. In the simulations of 
Table 6 Mg — 1.5 ■ 10^^ is the total mass in stars , with a small scatter between different 
runs. Thus SFR(/^b) ~ 15 ± Idh-^MQyr'^ / 10^'^ Mq for simulations I and III, while SFR(i^s) 
~ 66 ± 22h'~'^ MQyr~^ / 10^"^ Mq for simulation run II. These values can be compared with 
SFR(Coma) ^ 0.25h-'^MQyr-^/10^^MQ found by Kodama & Bower (2000) for the Coma 
cluster at 2; ~ 0.1. For the cosmology chosen by Kodama & Bower (fi^ = 0.2 and h — 0.5 
without a cosmological constant), at this redshift the Coma cluster has an age comparable 
to that of the simulated clusters at 2; = 0. Therefore simulations 1 and 111 have a cluster 
SFR which is consistent with the observed values. The simulation 11, with the KWH method 
and = 0.1, has instead an SFR clearly above the observed limits. For the high resolution 
runs, the cluster SF rates arc above the observed limits for Coma, with simulation IH being 
marginally consistent at a 2a level. The Coma cluster has been chosen for comparison 
because it has a measured cluster SFR and the estimated values for the mass components 
match those of ACDMOO. However the Coma cluster does not show a cooling flow activity 
and the cluster SF rates of the simulations would have been much smaller if ACDMOO did 
not have a cooling instability. 



4. Conclusions 

In this paper I have analyzed how the gas and X-ray properties of clusters of galaxies 
estimated from hydrodynamical SPH simulations are affected when radiative cooling is in- 
cluded. It has been found that in order to get reasonable results the inclusion of cooling 
cannot be decoupled from a prescription for converting cold gas particles into stars. The 
final results depend on the star formation prescription adopted and the numerical resolution. 
When the cooling gas particles arc converted into stars according to the NW prescription, 
the final cluster profiles are found to be remarkably stable under changes in the numerical 
resolution. 

This is achieved by using for the individual SPH cluster simulations a number of gas 
particles A^^^ ^ 20, 000 and a gas softening parameter Sg ^ 20h^^Kpc. Above this numerical 
threshold, estimates of the final cluster properties can differ among different runs by a factor 
^2. It must be stressed that these conclusions are vahd for the cluster sample studied 
here. As discussed in §2, numerical simulations of clusters less massive than those of Table 1 
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must take into account the metallicity dependence of the coohng function. In this case, 
the convergence of X-ray variables may require a numerical resolution higher than those 
considered here. 

If the KWH star formation prescription is adopted, with a star formation efficiency 
parameter = 1, the stability of the final results with respect to the numerical resolution 
of the simulation is satisfied for the same range of numerical parameters given above for the 
runs with the NW prescription. The KWH simulations with c^, = 0.1 have been found to give 
final results which are much more dependent on the simulation numerical resolution. For the 
KWH runs, the final differences in the gas density profiles at r ^ lOOKpc are a consequence 
of the different SF histories, which depend on c*. These differences strongly affect the X- 
ray luminosity Lx, which is then the simulation variable most sensitive to the value of 
c^. The results demonstrate that, for the same simulated cluster, different SF algorithms 
yield final gas distributions with differences locahzed at the cluster center. Global cluster 
properties, such as the total mass in stars, are robust to different SF prescriptions, while 
X-ray luminosities can differ by large factors. 

A relevant difference for the simulation runs with the NW prescriptions, with respect to 
previous simulations, is the flatness of the profiles. For the ACDMOO cluster the gas profile 

shows a core radius of Vc — 20kpc and the temperature profile is almost flat for r ^ lOOKpc. 
This is at variance with what expected for a cluster with a cooling flow, but the ACDMOO 
cluster has at its center a cooling time r,,(0) ^ 1/3 of the universe age tu ~ 13.5Gyr. 
The other two clusters which experienced a cooling flow (ACDM39 and CHDM39), have 
Tc{0) ~ IGyr and temperature profiles which decline by a factor ~ 2 within the ~ 200Kpc 
cluster central regions. 

A comparison of the temperature profiles with the results of Pearce et al. (2000) and YJS 
is problematic because these authors adopt a phcnomcnological prescription for removing 
cold gas particles from SPH estimates. In the case of YJS the temperature profiles include 
the cold gas population and show a steep dechne at the cluster centers. The profiles of the 
Cool-I-SF simulation of Lewis et al. (2000) can be compared with those of the corresponding 
KWH run with = 0.1 and do not show inconsistencies. Therefore the profiles obtained 
for the NW run are not inconsistent with previous simulations. The temperatiuc profiles 
show a decline between the cluster central regions and the virial radii, but the observational 
evidence for temperature gradients is controversial (Markevitch et al. 1998; Irwin Bregman 
& Evrard 1999; White 2000). 

To summarize, SPH hydrodynamical simulations of clusters of galaxies with radiative 
cooling and suitable star formation algorithms have been proved to be numerically stable, 
giving cluster X-ray properties which satisfy a set of observational constraints. An applica- 
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tion of the numerical schemes adopted to a simulated cluster sample can be used to reliably 
predict the evolution of the cluster X-ray luminosity and temperature function in different 
cosmological models. X-ray cluster surveys from the XMM mission can then provide strong 
constraints on the allowed cosmological background parameters, by comparing cluster data 
with simulation results. 
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Fig. 1. — Radial dependence of the gas density and temperature profile at z = in sim- 
ulations with radiative cooling. The left panel is for the test cluster ACDMOO, the right 
panel for ACDM39. In each panel, the upper plot is for densities and the lower plot is for 
gas temperatures. The density is in units of the critical density. The different curves are 
for integrations with different numbers of particles and different softening parameters (see 
Tables 2 & 3). The simulation run with index —00 is the integration without cooling. 
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Fig. 2.— As in Fig. 1, but for the test clusters CHDMOO and CHDM39. The numerical 
parameters of the simulations are given in Tables 4 and 5. 
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Fig. 3. — Final density and temperature profiles for the simulation runs including radiative 
coofing and star formation. The gas is allowed to convert into stars according to the NW 
star formation prescription (see text). The simulations are the same shown in Fig. 1, with 
the numerical parameters reported in Tables 3 & 5. The index of the simulations is clOO — k 
or cl39 — k, with k = j + 5 and j = 1,5 is the index of the cooling runs in Tables 2 
& 5. The simulations with index k = 11 have the same parameters as the k = 10 runs, 
but with the gravitational tolerance parameter 6=1 and quadrupole corrections enabled. 
The other simulations were performed with 9 = 0.7 without quadrupole corrections. In the 
density plots, the lines with a steeper slope (~ —3) show the density behavior of the star 
component. 
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Fig. 4. — As in Fig. 3, but for the cluster simulations shown in Fig, 2. 
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Fig. 5. — The cooling time Tc as a function of radial distance for the four test clusters. Tc 
is defined according to equation 6. In each panel, Tc is plotted for three different tests; the 
continuous line is the case with no cooling, the dotted line is the pure cooling test simulation 
with the highest resolution (—05), the short-dashed line is the equivalent simulation run but 
with gas particles being allowed to undergo star formation. 
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Fig. 6. — Mass accretion rates in the four test clusters for the same simulation tests shown 
in Fig. 5. In each radial bin, spherical averages for M(r) = Anpgr^Vr have been plotted only 
for negative values of the radial infall velocity Vr- 
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Fig. 7. — Final cluster X-ray luminosities as a function of Eg for the simulation runs performed 
for the four test clusters. In each panel Lx is shown for the simulations with gas cooling 
(open symbols) and for those also including star formation (filled symbols). For these test 
runs, values of the numerical inputs are reported in Tables 2, 3, 4 & 5. Simulation clusters 
with index clOO — j or cl39 — j correspond to the following symbols: diamond j = 00; circle 
j = 01,02; triangle j — 03,04; square j = 05. The filled symbols are the simulations with 
index j + 5. 
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Fig. 8. — Plots showing several cluster properties in simulation runs for the same cluster, but 
with different SF prescriptions. The test cluster is ACDMOO, and the simulation parameters 
are those of clOO — 11 (see Table 3). The different curves are for different SF methods or 
parameters, as reported in Table 6. The continuous line refers to the NW method, the 
others to KWH with different and 5^ ( c and e in the bottom right panel). Top left: star 
formation rate as a function of time. Bottom left: X-ray luminosity versus time. Top right: 
final radial density behavior for the gas and star components; for the sake of clarity the gas 
densities have been shifted downwards by a factor of 100. Bottom right: Radial temperature 
profiles. The dot-dashed line in the two left panels is for a simulation run with the NW 
method, but with an SN explosion energy of Esn — lO^^erg. 
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Fig. 9. — Baryonic fractions versus radius for the same simulations shown in Fig. 8. Top left: 
ratio of the total mass of gas within the radius r, to the cumulative cluster mass in units of 
the universal baryonic fraction. Top right: as in the left panel but for the ratio of (star+gas) 
mass to the total cluster mass. Bottom left: ratio of the star mass within radius r to the 
mass of gas within r. Bottom right: total mass produced in stars as a function of redshift. 



-43- 




0.1 0.2 0.3 0.4 0.5 

r [Mpc] 



Fig. 10. — Spherically averaged iron abundances Z — Fe/H as a, function of radius at the 
final epoch in units of the solar value 4.68 ■ 10~^. The different hues correspond to the 
simulations of Fig. 8. 
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Fig. 11. — For the first three runs I, II, III of Table 6, the simulation results of Fig. 8 are 
compared with the corresponding high-resolution runs IH, IIH and IIIH. The H simulations 
have the same SF parameters as the parent simulations, but the numerical parameters are 
given in the last row of Table 3. Left panels : the thick lines correspond to the high resolution 
simulations. Right panels: to facilitate a comparison with the high resolution results, the 
radial profiles of Fig. 8 have been shifted downwards by 10*^, A; = 2 for densities and k — 1 
for temperatures. 
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Table 1. Properties of the simulated clusters 



cluster 


M200 


^200 


CTl 


CHDMOO 


1.4 • 10^^ 


1.83 


1500 


CHDM39 


4.4 • 10" 


1.25 


1000 


ACDMOO 


1.2 • IQis 


1.98 


1200 


ACDM39 


4 ■ 10" 


1.37 


800 


Note. — 


Reference 


values at 


z = 



for the four clusters used in the numerical 
tests. M200; cluster mass within r2()() in 
units of H^^Mq, r2()o is in units of h^^ 
Mpc, (Ti is the central 1-D dark matter 
velocity dispersion in Kmsec~^ . M200 is 
defined as M200 = {47r/3)nriiPcAcr|()(), 
with Ac = 187Q^'^^ for a flat cosmol- 
ogy- 



Table 2. Numerical parameters of the simulations with cooling for 

the ACDMOO cluster 



ACDMOO 




m 


b 

9 


m 


c 
d 


N 

9 


^/ 


nJ 


93 


0'* 




clOO-00 


56 


3.01 ■ 


10"' 


2.64 


10" 


5503 


6295 


16463 


0.7 


F 


10. 


clOO-01 


56 


3.01 ■ 


10"' 


2.64 


10" 


5503 


6295 


16463 


0.7 


F 


10. 


clOO-02 


28 


3.01 ■ 


lOlO 


2.64 


10" 


5503 


6295 


16463 


0.7 


F 


10. 


clOO-03 


21 


1.45 • 


IQlO 


1.28 


10" 


11480 


14208 


35408 


0.7 


F 


10. 


clOO-04 


15.4 


1.45 • 


IQlO 


1.28 


10" 


11480 


14208 


35408 


0.7 


F 


10. 


cIOO-05 


21 


7.47 


• iqs 


6.57 


IQio 


22575 


25391 


67388 


0.7 


F 


19. 



Note. — Simulation parameters of the five test runs for the ACDMOO cluster. clOO — 00 is 
the reference case with no cooling, talcen from VGB. gravitational softening parameter for 
the gas in Kpc. mass of the gas particles in h~^MQ (the cosmology is for = 0.3 

and h = 0.7). : mass of the dark particles. number of gas particles inside the Lc/2 
sphere at 2 = Zi„. " : as in the previous column but for dark particles. ^ : total number of 
simulation particles, including those in the external shell of radius Lc- value of the treecode 
gravitational tolerance parameter. : gravitational quadrupole corrections F = disabled, T = 
enabled. ' : initial redshift for the simulation. 



Table 3. Numerical parameters for the simulations with cooling 
and star formation for the ACDMOO cluster 



ACDMOO 




m 


a 


ma 




Nd 


Nt 


e 


Q 




clOO-06 


56 


3.01 • 


IQlO 


2.64 


10" 


5503 


6295 


16463 


0.7 


F 


10. 


clOO-07 


28 


3.01 ■ 


IQio 


2.64 


10" 


5503 


6295 


16463 


0.7 


F 


10. 


clOO-08 


21 


1.45 ■ 


10"' 


1.28 


10" 


11480 


14208 


35408 


0.7 


F 


10. 


clOO-09 


15.4 


1.45 ■ 


10"' 


1.28 


10" 


11480 


14208 


35408 


0.7 


F 


10. 


clOO-10 


21 


7.47 


109 


6.57 


lOlO 


22575 


25391 


67388 


0.7 


F 


19. 


clOO-11 


21 


7.47 


10^ 


6.57 


lOio 


22575 


25391 


67388 


1.0 


T 


19. 


clOO-llH 


10.5 


2.45 


109 


2.12 


lOio 


69599 


74983 


204799 


1.0 


T 


29. 



Note. — As in Table 2, simulation parameters of the test runs including cooling and star 
formation for the ACDMOO cluster. The index of the run is clOO — k, with k = j + 5 and 

J = 1,5 is the index of the cooling simulation. The numerical parameters of the cooling and 
star formation runs are those of the corresponding cooling simulation with index j. The last row 
gives the numerical parameters for the high-resolution runs used to test different SF methods. 
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Table 4. CHDMOO 



CHDMOO 




mg 




a 
d 


^9 




Nt 




clOO-00 


50 


2.28 ■ lO^" 


3.57' 


10" 


5551 


13038 


27971 


4.8 


clOO-01 


50 


2.28 • lO^" 


3.57' 


10" 


5551 


13038 


27971 


4.8 


clOO-02 


25 


2.28 • IQl" 


3.57' 


■IQll 


5551 


13038 


27971 


4.8 


clOO-03 


15 


l.l-lOio 


1.73 


10" 


11507 


29038 


59093 


4.8 


clOO-04 


11 


l.l-lQio 


1.73 


10" 


11507 


29038 


59093 


4.8 


clOO-05 


15 


5.5 • 10^ 


8.9 ■ 


lOio 


22575 


50726 


112512 


9 



Note. — As in Table 2, but for CHDMOO. " : and Nj, are the total 
cold+hot values. 
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Table 5. ACDM39 and CHDM39 



ACDM39 CHDM39 





£9 


m 


a 




Eg 


m 


9 




C139-00 


42 


1.47 


lOio 


1.3-10" 


50 


7.5- 


109 


1.19 - 10" 


C139-01 


42 


1.47 


lOio 


1.3 ■ 10" 


50 


7.5 - 


109 


1.19 ■ 10" 


C139-02 


21 


1.47 


lOio 


1.3 ■ 10" 


25 


7.5 - 


109 


1.19 - 10" 


C139-03 


14 


7.2 ■ 


109 


6.37- 10^" 


15 


3.7 ■ 


109 


5.75 - 10^° 


C139-04 


10.5 


7.2- 


109 


6.37 - IQio 


11 


3.7- 


109 


5.75 - 10i° 


C139-05 


14 


3.7- 


109 


3.22 ■ IQio 


15 


1.89 


•109 


2.96 - IQi" 



Note. — Simulation parameters of the test runs for the clusters ACDM39 and 
CHDM39. The number of particles and initial redshifts are the same as for the 
00 clusters. 
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Table 6. Simulation runs with different star formations parameters 



method 


c* 


£* 




'T'a 
em 


rpb 

m 


r=^(< 50Kpc)/Tem 


run 


NW 


1 


1/2 


1 


8.3 


5.66 


0.08 


I 


KWH 


0.1 


1/2 


1 


1.9 


5.5 


0.8 


II 


KWH 


1 


1/2 


1 


8.96 


5.77 


0.2 


III 


KWH 


0.1 


1/3 


1 


2.05 


5.42 


0.87 


IV 


NW 


1 


1/2 


0.1 


7.88 


5.54 


0.08 


V 



Note. — SF parameters for the simulations used to compare difTerent models 
of SF. The test cluster is ACDMOO and the numerical parameters axe those of 
clOO — 11 (see Table 3). The different SF parameters are defined in §2, ssn is 
the SN explosion energy in units of 10^^ erg. " ; emission weighted temperature 
in KeV units at z = 0. : mass-weighted temperature at 2: = 0. irelative 
contribution to the total emission-weighted temperature from the r = 50Kpc 
cluster inner region. 



